Prom concentration profiles to polymer osmotic equations of state 
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We show that equilibrium monomer and centre-of-mass concentration profiles of lattice polymers 
in a gravitational field, computed by Monte-Carlo simulations, provide an accurate and efficient 
road to the osmotic equation-of-state of polymer solutions, via a straightforward application of the 
hydrostatic equilibrium condition. The method yields the full equation of state over a wide range 
of concentrations from a single simulation, and does not suffer from significant finite size effects. It 
has been applied to self-avoiding walk polymer chains with nearest neighbour monomer attractions, 
from the good solvent to the theta solvent regimes. The consistency of the method has been carefully 
checked by varying the strength of the gravitational field. 
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I. INTRODUCTION 

While the osmotic equation of state (e.o.s.) of off- 
lattice models of polymer solutions or melts can be 
readily computed from the contact theorem, according 
to which the osmotic pressure is proportional to the 
monomer density at a hard walli, the task is signifi- 
cantly more arduous for polymers on a lattice. For short 
chains or low polymer concentrations the chemical po- 
tential may be calculated from the insertion probability 
of a test chain- and the pressure then follows by stan- 
dard thermodynamic integration. The insertion method 
can be extended to solutions of longer polymers pro- 
vided configurational bias Monte-Carlo algorithms are 
used^a. An alternative method, the repulsive wall ther- 
modynamic integration (RWTI) method, which remains 
efficient at high polymer concentrations, or in the melt, 
was proposed by Dickman 5 , It extends the contact the- 
orem approach, but requires several simulation runs for 
increasing wall-monomer repulsion and subsequent ther- 
modynamic integration for each state point, rendering 
the method rather cumbersome. Moreover a recent anal- 
ysis by Stukan et. al& revealed that the RWTI method 
is prone to large finite-size effects, particularly at high 
polymer concentrations, which can only be overcome by 
switching to grand-canonical ensemble simulations. 

In this paper we show that the e.o.s. of dilute or semi- 
dilute polymer solutions is much more efficiently com- 
puted by subjecting the polymers to a strong gravita- 
tional field. The resulting sedimentation profile of the 
polymer solution then leads directly, via the hydrostatic 
equilibrium condition, to the osmotic e.o.s. over a wide 
range of concentrations from a single Monte-Carlo (MC) 
run. The method extends an idea which has been success- 
fully applied to concentrated dispersions of rigid colloidal 
particles (e.g. spherical^ or plate like&) to the case of 
flexible polymers; it applies to on- and off-lattice models 
alike. 



II. MONOMER AND CENTRE-OF-MASS 
CONCENTRATION PROFILES 

Sedimentation equilibrium of macromolecular solu- 
tions or colloidal dispersions arises from the balance be- 
tween gravity, which pulls particles to lower altitudes z, 
and the entropic drive toward homogeneity, and is char- 
acterised by a concentration profile p(z). If m is the 
buoyant mass of the particles, the sedimentation length 
is defined by £ = fc^T/mg, where g is the acceleration 
of gravity and k B T the thermal energy. For compact, 
micron-sized colloids, C is typically of the order of a 
few particle diameters, but for the much lighter fractal 
polymer chains, C is very large under normal gravity, so 
that the solution remains nearly homogeneous in prac- 
tice. However, in simulations g can be tuned to induce a 
measurable modulation of the concentration profile. 

We consider systems of N polymer chains, each of L 
monomers (or segments) living on a cubic lattice of spac- 
ing a. Each site hosts at most one monomer, corre- 
sponding to the self-avoiding-walk (SAW) model, while 
non-connected nearest neighbour segments interact with 
energy — e. The dimensionless inverse temperature is 
P* = 0e = e/k B T;/3* — corresponds to the ather- 
mal SAW limit, while the 0-solvent regime is reached for 
(3* « 0.265 for L=500 chains^. Each monomer experi- 
ences the gravitational energy —mgz, and the gravita- 
tional coupling constant is the dimensionless ratio A m — 
a/Cm — mga/krjT '. Sedimentation equilibrium is charac- 
terised by the monomer concentration profile p m (z), i.e. 
the mean number of monomers per unit volume at alti- 
tude z. A more coarse-grained representation focuses on 
the centre-of-mass (CM) of each polymer coil, of charac- 
teristic dimension equal to the radius of gyration R g . The 
gravitational field acts on the CM, and the corresponding 
sedimentation length is £ cm = k B T/Mg = Q m /L where 
M = Lm is the total mass of the polymer. The relevant 
gravitational coupling constant is 

_ R g _ MgR g R g 1+v 

Ccm K B 1 a 

where R g ~ aL u with v the Flory exponent (y ~ 0.59 in 



2 



good solvent and v = 0.5 in solvent). At sedimenta- 
tion equilibrium, the CM concentration profile is p cm {z), 
which satisfies the normalisation, 

/>OC 

/ p cm (z)dz = n s , (2) 
Jo 

where n s = N/A is the number of polymers per unit area 
of an x-y cross-section of the sedimentation column. 

Two limits of the concentration profile are known 
explicitly. First, the system of independent (free) 
monomers, i.e. L = 1 polymers, reduces in the SAW 
limit (j3* — 0) to a single occupancy lattice gas in a grav- 
itational field. The corresponding concentration profile 
/) m W is easily calculated to be: 



where p is the monomer chemical potential. At suffi- 
ciently high altitudes (z — > oo) the profile follows the 

barometric law for monomers: pm(z) ~ exp(— (3mgz) ~ 
exp(— z/Cm)- Upon introducing connectivity constraints, 
however, the concentration profile contracts significantly. 
At high altitudes z, where the polymer-polymer interac- 
tions are negligible, the concentration profile takes the 
form 

p m {z) oc p cm {z) - exp(~(3mLgz) ~ exp(-z/( cm ), (4) 

following the barometric law for polymers. In other 
words, it is contracted by a factor L over a system of un- 
connected monomers. This follows from the well-known 
1 /L reduction of the ideal component of the osmotic pres- 
sure of polymer solutions P td — p m /L, where p m is the 
bulk monomer concentration. In fact, if the e.o.s. and 
hence the free energy density of a homogeneous (g = 0) 
polymer solution is known as a function of bulk poly- 
mer concentration p and temperature, the full concen- 
tration profile in the weak modulation limit A <§; 1 can 
be easily calculated within the local density approxima- 
tion (LDA)i. Here we adopt the opposite point of view. 
We solve the inverse problem, and as explained in the 
next section, we extract the unknown e.o.s. from concen- 
tration profiles computed by MC simulations. 

The simulations were carried out for the cubic lat- 
tice model defined earlier, in a simulation box of dimen- 
sions l x x l y x l z in units of a. The square base of area 
A — l x x l y was periodically repeated in the x and y di- 
rections; most runs were for l x — l y = 100. The vertical 
dimension l z was chosen such that for reasonable val- 
ues of the gravitational coupling constant A cm , the poly- 
mer concentration at the highest altitude was negligible 
compared to the density at the bottom. Most simula- 
tions were carried out for N — 1600 polymers of length 
L = 500, using pivot, translation, reptation and con- 
figurational bias Monte-Carlo moves^. Starting from an 
initial homogeneous configuration, the system was equi- 
librated until the downward drift of the overall CM of 




FIG. 1: Monomer and centre of mass profiles for TV = 1600 
SAW polymer chains of length L = 500, shown for two val- 
ues of the gravitational coupling constant A C m. (a) Profiles 
plotted on a linear scale, (b) The same profiles, but now 
plotted on a logarithmic scale. The dashed lines denote the 
barometric law @, valid at low densities. 

the system stopped. The profiles p m {z) and p cm (z) were 
then calculated from altitude histograms averaged over 
several million configurations (the statistics are of course 
L times better for the monomer than the CM profiles.) 
For given N, L and /?*, simulations were carried out for 
several gravitational couplings A cm Local polymer con- 
centrations Pcm(z) are expressed relative to the bulk over- 
lap concentration p* = 3/4irR g 3 where R g is chosen to 
be the temperature dependent radius of gyration at zero 
concentration. For a given N, higher values of X cm are 
required to achieve higher polymer concentrations near 
the bottom of the simulation cell. 
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FIG. 2: Concentration profiles for N=1600 polymers of length 
L=500, at A cm =0.842, and /T=0. The dashed line is the 
convoluted profiles obtained from equation ©. 

Examples of concentration profiles p m (z) and p cm {z) 
for (3* — 0, a model for polymers in a good solvent, are 
shown in Figure ^ As expected, the profiles steepen 
when gravity, i.e. X cm increases. The CM profiles exhibit 
a marked first adsorption layer at the bottom, which is 
much smaller in p m (z) profilesii, but beyond that layer 
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the two profiles coincide within the statistical errors of 
p cm (z). The CM layering and preceding "correlation 
hole" in p cm (z) (which is more clearly apparent in Figure 
[2J may be traced back to the effective wall/ CM repulsion 
of entropic originAi. The logarithmic plots of the profiles 
reveal linear behaviour at high altitudes with slopes of 
— 1/Ccm, in agreement with the asymptotic barometric 
behaviour, thus providing a direct check on the conver- 
gence of the MC simulations. The close agreement be- 
tween p m {z) and pcm{z) beyond the CM adsorption layer 
is of course a consequence of the polymer connectivity. 
The agreement may be quantified by assuming that the 
monomer/CM form factor (which describes the distribu- 
tion of monomers around the CM in a polymer coil) is 
independent of the local polymer density and taken to 
be that appropriate for an ideal (Gaussian coil) polymer, 
namely: 



u C hi{r) = —e 



(~3r 2 /2R g 2 ) 



(5) 



y 2 + z 2 is the monomer/CM distance. 



where r — x 
The approximate relation between p m (z) and p cm {z) 
then simply given by the following convolution: 

/>oo 

p m {z) = / u>cm(\z - z'\)p cm (z')dz' . 



(6) 




FIG. 3: Theta solvent monomer and centre of mass concen- 
tration profiles. In the dilute regime the theta solvent e.o.s. 
resembles that of ideal polymers, and so follows the baromet- 
ric law for higher densities than polymers in good solvent do. 

An example is shown in Figure [3 The agreement be- 
tween the "exact" p m (z) and the approximation is seen 
to be good, except at distances less than R g from the wall 
where the internal structure of the polymers is expected 
to be distorted compared to its bulk behaviour. 

Concentration profiles p m {z) and p cm {z) under 9 sol- 
vent conditions(/3*=0.265) are shown in Figure[3]for two 
gravitational couplings X cm . The logarithmic plots show 
that the profiles reach the asymptotic barometric law ear- 
lier than in the SAW case, but begin to show important 
deviations from the asymptotic limiting law of Eq. Pfl 
at p p*. In other words polymers behave ideally in a 
^-solvent in the dilute regime p < p* onlyiS. 



III. FROM CONCENTRATION PROFILE TO 
OSMOTIC PRESSURE 

In order to extract the osmotic e.o.s. from the concen- 
tration profiles of section 2, we proceed as in reft For suf- 
ficiently slowly varying profiles, i.e. for sufficiently weak 
external field, the LDA holds and the Euler-Lagrange 
equation associated with the minimisation of the free en- 
ergy functional with respect to the concentration profile 
Pcm(z) leads back to the macroscopic equation of hydro- 
static equilibrium^: 



dP(z) 
dz 



-Mgp cm (z) 



(7) 



where P(z) is the local osmotic pressure at altitude z. 
Integration of Eq. (JTJ yields: 



i r 

(5P(z) = — / 

Scm J z 



p cm (z')dz' 



(8) 



Thus P(z) and p cm (z) are known as functions of altitude, 
with elimination of z leading to the desired e.o.s. of the 
bulk polymer solution P = P(p cm ,T). Note that the 
LDA approximation is not expected to be accurate near 
the hard wall where p cm (z) varies rapidly within the first 
adsorption layer, so that the lower integration limit in 
Eq. JHl should be taken at z > R g , i.e. beyond the peak 
of Pcm(z), where p cm (z) and p m {z) become indistinguish- 
able. 

Eq. (JHl has an obvious intuitive interpretation: the sys- 
tem relaxes to a density p(z) such that the local pressure 
P(p(z)) counteracts the weight of the polymers above. 
Conversely, when the e.o.s. is known, Eq. Q can be di- 
rectly integrated to calculate the concentration profile 7 . 
For example, in the semi-dilute regime, where the os- 
motic equation of state obeys the well known scaling 
law (3P(p)/p ~ 7(p/p*) 1 /( 3,y ~ 1 Ui, with 7 a dimensionless 
constant that depends on polymer details, the resulting 
concentration profile is given by: 



pr*{z) = po 1- ^ 



3//-1 



37J/ Rg 



(9) 



where po is the maximum polymer density, reached at 
z = 0, and follows from inverting the equation 



P(po) = Mgn s 



(10) 



if the wall-induced layering at the bottom of the con- 
tainer is ignored. Eq. @ provides an accurate fit to the 
density profile in the semi-dilute regime, but at higher 
altitude the analytic profile p(z) fails to cross over to the 
barometric law as illustrated in Figure (JJJ). 

We have carried out the above inversion procedure for 
/3*=0, 0.1, 0.2, and 0.265 (9 solvent). The resulting e.o.s. 
P(p, T) should be independent of the gravitational cou- 
pling used in the MC simulations provided X cm is not too 
large, i.e. the gravitational field is not too strong. This 
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FIG. 4: Concentration profiles, compared to the analytic form 
from equation ©. Note the excellent agreement in the semi- 
dilute regime and the expected disagreement in the dilute 
regime where the barometric law takes over. 

was checked explicitly by carrying out the inversion pro- 
cedure for two different values of A cm ; the resulting e.o.s. 
always turned out to be practically identical, with very 
small differences (typically less 1% over the whole con- 
centration range) providing an estimate for the sum of 
systematic and statistical errors. However, if the applied 
gravitational field is too strong (A cm > I) the result- 
ing profiles vary too rapidly with z for LDA to remain 
accurate. This reflects itself in the failure of the result- 
ing compressibility factor Z=/3P/p to go over to its ideal 
gas limit as p — > 0, as illustrated in the inset of Fig- 
ure Our results for Z(p, T) are plotted in Figure [S] 
along four isotherms, and compared to the predictions of 
MC simulations of bulk polymer solutions 12 based on the 
RWTI procedure^. The agreement is seen to be excellent 
throughout, tn fact some apparent discrepancies between 
our earlier datai^, based on the RWTI method and those 
obtained with the present method were resolved when 
contentious RWTI high concentration points were re-run 
with considerably improved statistics. 

In an effort to detect possible finite size effects, we 
repeated some of the simulations with a four times larger 
base area. Any observed differences in the resulting data 
were within statistical noise. 

A related consistency test of the LDA-based inversion 
procedure is provided by the scale invariance property 
of the LDA concentration profile, according to which, 
along any isotherm, the re-scaled profiles p C m(z/£ cm ) 
or p m (z/Ccm) depend only on the dimensionless prod- 
uct n s X cm a 2 . This means that as long as LDA applies, 
the re-scaled profiles should be identical if the number 
of polymers is divided by a given factor, provided the 
gravitational coupling A c is multiplied by the same fac- 
tor. An illustration of this scale invariance is shown in 
Figure HO The equation-of-state extracted from the den- 
sity profiles are in perfect agreement. Deviations from 
scale invariance would signal the breakdown of the LDA 
assumption, providing a diagnostic for the consistency of 
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FIG. 5: Equation of state isotherms calculated by inver- 
sion of concentration profiles (lines) compared to earlier 
simulations 12 using the RWTI method (symbols). Whereas 
the former need about 50 independent simulations per 
isotherm, the present method requires only one simulation 
per isotherm. The solid lines are for runs \ cm < 1, the dotted 
(3* = line is at A cm = 4.21 and the dashed /3* = 0.1 line is 
at A cm = 1.97. 




FIG. 6: (a) Three different concentration profiles from simu- 
lations with different numbers of polymers and gravitational 
strengths, but with the product n s X cm a 2 held constant. Note 
that the equations of state (b) generated from the concentra- 
tion profiles are practically identical, as are the concentration 
profiles (c) when scaled by the gravitational length. 



IV. EFFICIENCY OF INVERSION METHOD 

By inverting the concentration profile in an imposed 
external field, we measure the full equation of state 
in a single simulation. This contrasts with the RWTI 
method 5 , where for each value of p a number of separate 
simulations are needed to perform thermodynamic inte- 
gration. For example, when we calculated the e.o.s. of 
L = 500 polymers as a function of solvent quality in 12 , 
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using the RWTI method, we needed about 5 simulations 
for each of the 10 statepoints along each isotherm. In the 
current work, we used only 1 simulation of similar size to 
each of the 50 simulations per isotherm carried out ini£. 
As demonstrated in Fig. |5J the results are of comparable 
accuracy, but achieved with considerably less CPU time. 

For a given A cm and number of polymers N, the max- 
imum density achieved under an external field can be es- 
timated from Eq. (jTU|l . If the e.o.s. scales as (3P(p)Rg ~ 
(p/p*) Q , then the maximum density pg scales as 

Larger po can be achieved by increasing A cm , although 
this is constrained by the LDA criterion, or by increas- 
ing N. For semi-dilute polymers a ss 2.3 in good solvent 
and a w 3 in a theta solven^^. To double po, the 
number of polymers N must be increased by a factor five 
for the former and eight for the latter solvent quality. 
On the other hand, for the RWTI method, doubling po 
simply means doubling N, irrespective of a. Of course 
to calculate an isotherm, this may also imply doubling 
the number of different p at which independent simula- 
tions must be performed. Furthermore, as pointed out by 
Stukan et. al. e finite size effects become more severe for 
larger polymer concentrations. As an example, they cal- 
culated the pressure of bond-fluctuation model polymers 
of length L = 20, using a box with a width l x = l y = 20 
and varying the distance l z between the two hard walls 
used for thermodynamic integration. At the rather high 
monomer packing fraction <fi — NL/V = 0.5 they used, 
the influence of the two hard walls was noticeable even up 
to distances of l z = 160, where the RWTI method overes- 
timated the pressure by about 3%. We performed compa- 
rable simulations for L = 20 SAW polymers at the same 
melt-like density. In a gravitational field we were able to 
reproduce the l z — > oo results of the RWTI method, but 
with the number of polymers N which the latter method 



would need for a size l z — 20 box only. Thus the less 
favourable scaling of our external field method with N is 
partially offset by less severe finite size effects, allowing 
significantly fewer particles to be used. 

V. CONCLUSIONS 

We have shown that the hydrostatic equilibrium 
method, which allows the e.o.s. of polymer solutions to 
be computed along an isotherm in a single simulation, 
which determines the polymer concentration profiles in 
a gravitational field, is both accurate and efficient. It 
provides the osmotic pressure as a function of concentra- 
tion with a computational effort which is a small fraction 
(typically under 5%) of that required using the RWTI 
method, since the latter requires a series of MC simula- 
tions at different densities, and several independent sim- 
ulations at each density are needed for thermodynamic 
integration. Moreover the finite size problems which can 
affect the RWTI method are insignificant in the present 
method, mainly because the simulated system is "open" 
at high altitudes, i.e. essentially extends to infinity in the 
z-direction. The hydrostatic equilibrium method applies 
equally well to on and off-lattice models of interacting 
polymers. The only apparent limitation of the method 
is that a large gravitational field, i.e. large values of A cm 
are required to achieve large densities. For too large X cm 
the underlying LDA becomes inaccurate, and the inver- 
sion procedure can lead to erroneous results. Hence the 
method is expected to be well adapted to dilute and semi- 
dilute solutions, while to reach polymer densities typical 
of polymers melts, several runs with increasing values of 
A cm are needed to cover successive and partially overlap- 
ping ranges of p/p*. 

We plan to extend the present method to determine 
the osmotic e.o.s. of more complex polymeric systems, 
including block copolymer solutions and melts. 

CIA would like to thank the EPSRC for funding, AA 
is grateful for the support of the Royal Society. 



1 J.K. Percus, J. Stat. Phys. 1976, 15, 423. 

2 A. Bellemans, E. Devos, J. Polym. Set. Symp. 1973, 42, 
1195, R. Dickman and C.K. Hall, J. Chem. Phys. 1986, 
85, 3023. 

3 J.I.Siepmann, I.R. McDonald, D. Frenkel, J. Phys. Con- 
dens. Matter 1992, 4, 679. 

4 D. Frenkel, B. Smit, "Understanding Molecular Simula- 
tion", 2 nd Ed. Academic Press, New York, 1996. 

5 R. Dickman, J. Chem. Phys. 1987, 87, 2246. 

6 M.R. Stukan, V.A. Ivanov, M. Miiller, W. Paul, K. Binder, 
J. Chem. Phys. 2002, 117, 9934. 

7 T. Biben, JP. Hansen, J.L. Barrat, J. Chem. Phys. 1993, 
98, 7330. 



8 M. Dijkstra, JP. Hansen, P.A. Madden, Phys. Rev. E. 
1997, 55 3044. 

9 R. Piazza, T. Bellini, V. Degiorgio, Phys. Rev. Let. 1993, 
71 4267. 

P. Grassberger, R. Hegger, J. Chem. Phys. 1995, 102, 
6881. 

1 P.G. Bolhuis, A.A. Louis, J. P. Hansen, E.J. Meijer, J. 
Chem. Phys. 2001, 114, 4296. 

2 C.I. Addison, A.A. Louis, J. P. Hansen, J. Chem. Phys. 
2004, 121, 612. 

3 R.H. Corey, M. Rubenstein, Polymer Physics (Oxford Uni- 
versity Press, Oxford, 2003 



